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1.  Introduction 

Spectral  methods  have  been  used  to  great  advantage  in  hydrodynamic  stability  calculations; 
the  concepts  are  described  in  Orszag’s  seminal  application  of  the  Chebyshev  tau  method  to 
the  Orr-Sommerfeld  equation  for  plane  Poiseuille  flow  in  1971  [1  ].  Orszag  discusses  both 
the  Chebyshev  Galerkin  and  the  Chebyshev  tau  methods,  but  presents  results  for  the  tau 
method,  which  is  easier  to  implement  than  the  Galerkin  method.  The  tau  method  has  the 
disadvantage  that  two  unstable  eigenvalues  are  produced  that  are  artifacts  of  the  discretiza- 
tion. The  occurrence  of  spurious  eigenvalues  has  been  discussed  by  several  authors,  c.f. 
[2, 3, 4, 5, 6, 7], 

In  this  note  we  present  an  extremely  simple  modification  to  the  Chebyshev  tau  method 
which  eliminates  the  spurious  eigenvalues.  We  first  study  a simplified  model  of  the  Orr- 
Sommerfeld  equation  discussed  by  Gottlieb  and  Orszag  [2].  We  consider  the  Chebyshev  tau 
method,  which  has  two  spurious  eigenvalues,  and  then  describe  a modification  which  elim- 
inates them.  Our  modification  is  motivated  by  considering  two  other  discretizations  of  the 
model  problem  which  also  have  no  spurious  modes:  a vorticity-streamfunction  reformulation 
of  the  Chebyshev  tau  method,  and  the  Chebyshev  Galerkin  method.  For  the  model  problem 
we  show  that  the  latter  approaches  are  equivalent,  and  that  both  reduce  to  our  modification 
of  the  tau  method.  We  also  remark  on  the  modified  Galerkin  method  formulated  by  Zebib 
[4,6].  Finally,  we  consider  results  for  the  Orr-Sommerfeld  equation,  where  our  modified  tau 
method  also  eliminates  the  spurious  eigenvalues.  The  simplicity  of  the  modification  makes 
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it  a convenient  alternative  to  other  approaches  to  the  problem. 


2.  The  Model  Problem 


In  their  1977  monograph  on  spectral  methods,  Gottlieb  and  Orszag  [2,  p.  143]  consider  the 
problem 


^zzt  = ^ 


(-1  < * < 1,  t > 0),  (1) 


with  the  boundary  conditions  $(±1,  t ) = $,(±1,  t)  — 0,  as  a simple  model  of  incompressible 
fluid  dynamics.  A normal  mode  of  the  form  T(2,  t ) = ip{z)  exp  at  gives  rise  to  the  eigenvalue 
problem 

*l>zzzz  = (-1<2<1),  (2a) 


with  the  boundary  conditions 


V>(±1)  = tG(±l)  = 0 


(2b) 


The  latter  problem  models  some  features  of  the  Orr-Sommerfeld  equation  [1],  and  can  be 
solved  exactly.  The  eigenvalues  of  (2)  are  given  by  the  values  a — — /r2,  where  either  f.i  = nir 
for  n = 1,2, . . .,  corresponding  to  the  even  eigenfunctions  with  x^>zz  = cos  fiz,  or  ^ is  a positive 
root  of  the  transcendental  equation  tan  p = corresponding  to  the  odd  eigenfunctions  with 
xjjzz  — sin  [iz. 


2.1.  The  Tau  Method 

Gottlieb  and  Orszag  [2]  show  that  a straightforward  application  of  the  Chebyshev  tau  method 
to  Eq.  (2)  gives  rise  to  positive  eigenvalues  whose  magnitudes  increase  rapidly  as  more  terms 
are  included  in  the  expansion.  Such  spurious  modes  may  be  discarded  by  inspection  for 
the  system  given  by  (2),  but  they  cause  severe  numerical  instability  in  the  time-dependent 
system  (1).  The  equations  for  the  Chebyshev  tau  method  are  obtained  by  writing 

N 

</>  = £ a„T„(z)  (3) 

n—0 
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where  Tn(z)  is  the  n-th  degree  Chebyshev  polynomial,  which  satisfies  Tn(cos9)  = cos  nO.  We 


then  have  [l] 

N N 


^zz  = a{n]Tn{z 

), 

(4a) 

71=0 

77=0 

where 

N 

Cna(n]  = 

E 

p(p2  - n2)apl 

(4b) 

p=7i+2 
P+tz  even 

cna^4'  = — 
24 

E pIp2(p2 

_4)2_ 

3n2p4  + 3n4p2  — n2(n2  — 4)2]ap, 

(4c) 

p=n+4 
p+n  even 

and  Cq  = 

2 and  cn  = 1 

(2) 

for  n > 0.  Note  that  these  expressions  give  0 = aN_ l 

- a(2)  - 

— aN  — 

a[N_3  — a;v_ 2 = aW-\  — aN  ■>  since  V’z*  and  V’zzzz  are  polynomials  of  degree  N — 2 and  N — 4, 
respectively.  Equations  for  the  coefficients  an  are  defined  by  the  relations 


frr  _ /rr 

dz 4 } “ <7(r?’  ’ 


i — 0, . . . , A — 4, 


where  in  the  Chebyshev  inner  product  we  have  (Tj,T^)  = (tt  !2)cj6jk.  The  resulting  tau 
equations  are  then 


<7  a 


(2) 


n = 0, . . . , N — 4,  (5a) 


with  the  four  boundary  conditions  that  follow  from  (2b), 

N N 

'£(±l)nan  = £(±l)”n2a„  =0.  (5b) 

n— 0 n=0 

The  equations  can  be  written  as  a generalized  eigenvalue  problem  Ax  — aBx , where  ,4  and 
B are  (N  + 1)  x (N  + 1)  matrices  whose  first  N — 3 rows  are  defined  by  Eq.  (5a).  The  last 
four  rows  of  A are  given  by  Eq.  (5b),  and  the  last  four  rows  of  B vanish. 
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We  illustrate  the  matrix  equations  schematically  for  the  case  N = 8: 
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here  the  non-zero  entries  are  denoted  by  the  symbol  V. 

This  system  may  be  solved  numerically  using  software  such  as  the  routine  RGG  from  the 
EISPACK  library  [8].  In  Table  I we  give  numerical  results  for  the  system  (5)  obtained  with 
RGG.  The  table  shows  the  values  of  the  first  and  fifth  negative  eigenvalues  <7i  and  and 
the  larger  of  two  positive  spurious  eigenvalues  <rmax.  These  results  reproduce  those  given  in 
Table  13.1  of  [2].  (Note  that  there  is  a typographical  error  in  [2]  in  the  results  for  <75  for 
.large  N.)  The  numbers  were  computed  in  double  precision  on  the  CDC  Cyber  205  at  the 
National  Institute  of  Standards  and  Technology  to  reduce  the  round-off  error.  For  single 
precision  calculations  on  the  Cyber  205  the  effect  of  round-off  becomes  noticeable  in  the  last 
two  digits  of  <7i  for  N > 25;  the  given  values  for  cr5  and  crmax  remain  the  same. 

Subroutine  RGG  provides  solutions  to  the  problem  written  in  the  form  (3 Ax  = aBx  [9],  so 
that  for  (3j  ^ 0 the  j-th  eigenvalue  is  given  by  crj  — ctj/ flj.  Solutions  with  f3  = 0 correspond  to 
the  null  space  for  the  related  problem  Bx  = fiAx,  where  n — 1/cr,  and  may  be  associated  with 
infinite  values  for  <7.  Eigenvalues  with  /3  = 0 occur  for  the  above  system  due  to  the  boundary 
conditions  (5b).  For  a given  value  of  N > 6,  four  eigenvectors  with  (3  = 0 are  obtained,  and 
the  other  N — 3 eigenvectors  have  (3^0.  Eigenvalues  with  (3  — 0 may  be  avoided  by  using 
the  four  equations  (5b)  representing  the  boundary  conditions  to  eliminate  the  variables  ayv- 3, 
ayv—  2,  ayv- 1,  and  ayv  from  the  system  as  described  in  [7].  This  eliminates  the  eigenvalues  with 
(3  = 0;  the  values  with  (3^0  are  unchanged.  The  numbers  aq  and  (3j  can  also  be  used  [10] 
to  give  an  estimate  of  the  sensitivity  of  the  computed  eigenvalue  to  perturbations  in  A and 
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B:  <Tj  is  ill-conditioned  if  the  condition  number  C(aj)  = 1/ \J\aj\1  + \(3j\2  is  large,  in  which 
case  crj  will  be  difficult  to  determine  using  low  precision  arithmetic. 

2.2.  A Modified  Tau  Method 

A useful  modification  of  the  tau  method  (5)  is  obtained  by  considering  instead  the  equations 


a(4)  = aa(2) 

u n u un  i 


n = 0, . . . , N — 4, 


i a 


with  the  four  boundary  conditions 

N 


N 


^(±1  )nan  — ]P(±l)nn2an  = 0, 


7b) 


n—  0 


n=  0 


where  we  have  defined 

Cnfli2’  = p(p2  - n2)civ.  (8) 

p=n  + 2 
p+Ti  even 

Note  that  a ^ is  obtained  from  a ^ by  simply  dropping  the  last  two  terms.  Thus,  the 
modified  matrix  eigenvalue  problem  is  obtained  from  the  tau  eigenvalue  problem  Ax  — aBx 
by  setting  the  last  two  columns  of  B to  zero. 

Numerical  results  for  the  system  (7)  are  also  presented  in  Table  I.  The  two  spurious  eigen- 
values that  are  obtained  using  the  standard  tau  method  are  eliminated,  and  the  remaining 
eigenvalues  are  computed  with  essentially  no  loss  in  accuracy.  The  numbers  were  computed 
in  double  precision  on  the  CDC  Cyber  205  at  NIST. 

The  motivation  for  our  modification  of  the  usual  tau  method  is  given  in  the  following 
two  sections. 


2.3.  Vorticity-Streamfunction  Formulation 

Gottlieb  and  Orszag  [2]  show  that  the  spurious  roots  in  the  Chebyshev  tau  method  are 
eliminated  if  the  function  ( = tpzz  is  introduced;  £ plays  the  role  of  vorticity  in  this  simple 
model.  In  this  case  the  model  equations  take  the  form 


(9a) 

tyzz  = C, 

(9b) 
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with 


VK±i) 

= V>*(±i)  = o. 

(9c) 

We  now  have  two  coupled  second  order  equations  for  ( and 

ip;  note, 

however, 

that  there  are 

four  boundary  conditions  on  ip  and  none 

on  (. 

If,  in  addition  to  Eq.  (3),  we  write 

N 

c = 

77=0 

(10) 

then  the  tau  equations  for  (9)  become 

e = 

&bn , 

n = 0, . 

1 

to 

(lla) 

a<2)  = 

' bm 

n = 0, . . 

N - 2 

..,1V  w, 

(lib) 

together  with  the  four  boundary  conditions 

N 

N 

Y2(±l)nan 

= £(±1  )"n2a„  = 

= 0. 

(11c) 

71=0  71=0 


These  equations  also  take  the  form  of  a generalized  eigenvalue  problem,  and  may  be 
solved  using  subroutine  RGG.  In  this  case  it  is  found  that  the  discretized  Poisson  equation 
(lib)  contributes  N + 1 eigenvectors  with  (3  = 0.  There  are  an  additional  four  eigenvectors 
with  (3  = 0 arising  from  the  terms  representing  the  boundary  conditions,  leaving  N — 3 
eigenvectors  with  (3^0.  For  this  system  all  eigenvalues  with  (3^0  are  negative;  the 
spurious  eigenvalues  have  been  eliminated. 

As  pointed  out  by  Gardner  et  al.  [7],  the  size  of  the  system  may  be  reduced  by  using  the 
last  two  equations  in  (11a)  to  eliminate  the  unknowns  i and  6jv,  and  then  using  (lib) 
to  eliminate  bj  for  j = 0,...,iV  — 2.  The  first  N — 3 equations  in  (11a),  together  with  the 
boundary  conditions  (11c),  then  constitute  a generalized  eigenvalue  problem  analogous  to 
the  system  (5),  but  without  spurious  eigenvalues  that  have  Re((j)  > 0.  The  eigenvalues  with 
(3  0 are  unchanged  by  the  elimination.  We  next  show  explicitly  the  difference  between 

formulations  (5)  and  (11).  If  we  introduce  the  general  notation 

N 

= j = 0,...,N, 

k= 0 
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where  the  elements  Y jk  follow  from  the  explicit  form  of  Eq.  (4b),  then  the  last  two  equations 
in  (11a)  give 

vbN_2  = b$  = r tv— 2,7V^7V}  crbw-z  — bj^_  x = T/v-s.tv-i&tv-i • 


Putting  these  expressions  for  6yv  and  6yy- 1 into  the  first  N — 3 equations  in  (11a),  and  then 
using  (lib)  to  eliminate  the  remaining  6;,  we  obtain 


'V'*'  r (2)  ( (2)  ^jV-l  (2)  Yj  N (2) 

2_,  I flj  - p %-3  - p aN-2 

k=  0 V 1 W-aTV-l  i 


TV— 3,7V— 1 i 7V-2.7V 

Since  u'v-i  — aN  = 0,  the  left  hand  side  may  be  written 
TV— 2 TV 


j = 0,...,lV  — 4. 


.(2) 


£ ?A2)  = £ i> <2>  = 


i = o jv. 


0 


fc=0 


Finally,  since  a|v-3  — ryv-3,TV-i<2TV-i  and  a^_2  = Yjv_2,NaNi  we  conclude  that  the  tau 
equations  (11a)  and  (lib)  are  equivalent  to 


<44)  = cra^]  - crYri'N-iaN-i  - a Y 


n,N^  N 


n 


= 0, . . . , iV  — 4. 


' 12' 


Since 


ai,2>  = 


TV 

^ ( rnfca* 

yfc=0 


the  resulting  system  is  precisely  equivalent  to  our  modified  equations  (7). 


2.4.  Chebyshev  Galerkin  Method 

The  Chebyshev  Galerkin  discretization  of  the  model  problem  (2)  is  also  free  from  spurious 
eigenvalues  [2].  In  this  representation,  we  write 

TV 

an<i>n{z),  (13) 

n=  4 

where  the  functions  <f>n  are  chosen  to  satisfy  the  boundary  conditions  <j>n  = d4>n/dz  = 0 at 
2 = ±1,  i.e. , 


4>n{z)  — Tn(z)  + ')noTo(z)  + ')n\Tx(z)  + 7n2^2(2)  + ln?,T2{z)  Tl  — 4,  . . . , N] 
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here 


7n0  — 


7n2 


\{n2  — 4),  72  even 


0, 


72  even 


0, 


odd 


7nl  — 


|(  n2  — 9),  72  odd 


— i-722,  72  even 


7n3  — 


0,  n even 

— Un2  — 1),  n odd. 


0,  72  odd 

The  Galerkin  equations  for  the  coefficients  are  defined  by  the  relations 

d4ip  N / ( d2ip 

dT2 


( ’ ~dz*  ) = 


j = 4,...,Ah 


(14) 


If  we  introduce  four  additional  coefficients  a0,  a*,  a2,  and  a3,  we  may  write  as  before 


N 


V>( 2 ) = n anTn(z), 


(15) 


n=0 


where 


It  follows  that 


N 


— 'y  ^ 0‘n'lnj i J — 0,  . . . , 3. 


71=  4 


d4lp  _d2Xp  rj-1  / ^\ 

.4  G J 2 ~ / ; en-^  n(^  b 


(16) 


n=  0 


d.?4 

where  en  = a^4*  — a a Gfi  The  Galerkin  equations  (14)  are  therefore  equivalent  to  the  iV  — 3 
equations, 

0 — G'  + 27joCo  + 7jiei  + lj2^2  + 7j363  j = 4, . . . , N,  (17a) 


plus  the  boundary  conditions 

N N 

^(±l)nan  = ^(±l)n72  2an  = 0 (17b) 

n= 0 n= 0 

which  hold  by  the  construction  of  a0,  a1?  a2,  and  a3. 

We  next  show  that  equations  (17a)  are  equivalent  to  (12).  Since  t p(z)  is  a polynomial 
of  degree  N,  it  follows  from  (16)  that  e^r  — e^-i  — 0,  e^_2  = — ctT/v^a^/v  and  e;v_3  = 
— crr/v-3,iv-iaA7-i-  The  last  four  equations  in  (17a)  may  then  be  used  to  obtain  e0,  e1?  e2 
and  e3  in  terms  of  and  cln-\\  a calculation  gives 


tk  — — crT k,N-iaN-i  — vFk'Na-N,  k = 0, ...  ,3. 

A further  calculation  then  shows  that 


27j0e0  + 7jlel  + 1]2&2  + 7j363  — ^ j,N-laN-l  + CT  T j'NaN 


j = 4, . . . , TV  - 4. 


Combining  these  last  two  expressions  with  (17a)  we  obtain 

en  = — oTnijv-iaiV-i  — cT n,NaN  n — 0, . . . , N — 4, 

which  reduce  to  (12). 

We  conclude  that  each  of  the  discrete  eigenvalue  problems  (7),  (11),  and  (17)  are  equiv- 
alent in  the  sense  that  eigenvalues  a with  f3  / 0 that  are  produced  by  each  formulation  are 
identical. 


2.5.  Other  Galerkin  Procedures 

Other  variants  of  the  Galerkin  formulation  are  possible.  For  example,  Zebib  [4,6]  introduces 
a different  basis  by  writing 

N-A 

(18) 


N-A 

hnUn{z), 

n—  0 


where  the  un(z)  are  linearly  independent  polynomials  of  degree  at  most  iV,  uniquely  char- 
acterized by 

dS 


, =Tn(z), 

Qj4, 

with  the  boundary  conditions  un  = dun/dz  = 0 at  z = ±1.  Zebib’s  original  approach  [4] 
produced  results  that  included  spurious  eigenvalues,  which  he  was  later  able  to  eliminate  [6]. 

The  functions  un(z)  for  n — 0, . . . , N — 4,  and  the  functions  (j>k  for  k = 4, . . . , N used 
in  the  previous  section  span  the  same  space,  i.e.  the  set  of  N- th  degree  polynomials  which, 
together  with  their  derivatives,  vanish  at  z = ±1.  It  follows  that  the  function  ip(z)  Eq-  (18) 
can  be  re-expressed  as  a linear  combination  of  Tn(z),  n = 0,...  ,iV,  as  in  Eq.  (15),  where 
the  coefficients  a3  are  linearly  related  to  the  constants  hn  above  and  satisfy  (17b).  Eq.  (16) 
holds  as  well. 

Using  Zebib’s  original  approach  [4]  on  the  model  problem,  we  take  the  inner  product 
of  the  equations  (16)  with  Tn  for  n = 0,...,jV  — 4 and  obtain  equations  which,  in  terms 
of  the  variables  are  precisely  equivalent  to  the  Chebyshev  tau  equations  en  = 0 for 
n = 0, . . . , N — 4. 

Using  Zebib’s  later  approach  [6]  we  take  the  inner  product  of  the  equations  (16)  with  un 
for  n = 0, . . . , N — 4.  Since  each  basis  function  un  can  be  expressed  as  a linear  combination 
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of  the  functions  0^,  we  obtain  a linear  combination  of  the  equations  (17).  The  corresponding 
matrix  equations  are  now  related  by  pre-  and  post-multiplication  by  invertible  matrices,  and 
the  spectrums  are  identical.  Thus  Zebib's  second  approach  produces  eigenvalues  which  are 
in  principle  identical  to  those  obtained  using  the  regular  Galerkin  formulation  outlined  in 
the  previous  section.  In  practice,  the  condition  numbers  for  the  eigenvalues  obtained  using 
Zebib’s  later  approach  are  large  compared  to  those  obtained  using  the  standard  Galerkin 
method.  Roughly  speaking,  this  may  be  attributed  to  the  fact  that  the  basis  un  is  not  as 
well-conditioned  as  the  basis  0^;  the  0*.  are  more  nearly  orthogonal  (cf.  [11]). 

These  results  generalize,  in  the  sense  that  expanding  in  any  basis  for  polynomials  of 
degree  N satisfying  the  boundary  conditions  still  produces  a function  which  satisfies  Eqs. 
(15)— ( 16).  Taking  the  inner  product  of  the  equations  with  Tn  will  produce  the  spectrum 
from  the  Chebyshev  tau  method,  with  spurious  eigenvalues.  Taking  the  inner  product  of  the 
equations  with  the  basis  functions  will  produce  the  spectrum  from  the  standard  Galerkin 
method,  with  no  spurious  eigenvalues. 


3.  Other  Model  Problems 


In  Table  II  we  list  eight  model  problems  of  second  through  sixth  order  which  were  discretized 
and  solved  using  the  Chebyshev  tau  method  with  N — 20.  Spurious  eigenvalues  were 
observed  in  three  of  these  cases.  The  first  was  considered  in  Section  2.  The  spurious  modes 
in  the  remaining  two  problems,  both  sixth  order,  were  eliminated  using  a modified  tau 
method  similar  to  that  proposed  in  Section  2.2. 

In  problem  7 we  use  the  modified  discretization 

a If*  = era ^ n = 0, . . . , N — 6, 


where 


cna(n' 1 = X]  p{p2-n2)ap- 

p=n+ 2 
p+n  even 


that  is,  the  last  four  columns  of  the  matrix  representing  are  set  to  zero.  In  problem  8 
we  use  the  discretization 


a 


(6) 

n 


a a L4' 


n = 0, . . . , N — 6 
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where 


cna'n‘>  = 77  S P[P2(P2  ~ 4)2  - 3 n2p4  + 3 n4p2  - n2(n2  - 4)2]ap; 
24 


d4)  = _ 


p=n+4 

p+n  even 


that  is,  the  last  two  columns  in  the  matrix  representing  a [T  have  been  set  to  zero.  In  each  of 
these  cases  no  significant  change  in  the  accuracy  of  the  computed  eigenvalues  was  observed. 


4.  Orr-Sommerfeld  Equation 

The  Orr-Sommerfeld  equation  for  plane  Poiseuille  flow  may  be  written  in  the  form  [1] 

[^zzzz  - 2 a2^zz  + a4^]/(iaR)  - ( U(z ) - s)(ipzz  - a2ip)  + Uzzip  = 0 

for  — 1 < z < 1,  where  U(z)  — 1 — z2  is  the  base  velocity,  a is  the  wavenumber  of  the 
disturbance,  R is  the  Reynolds  number  of  the  flow,  and  s is  the  temporal  eigenvalue.  The 
boundary  conditions  at  2 = ±1  are  ifj  = tpz  = 0.  For  a given  Reynolds  number,  the  flow 
is  stable  if  Im(s)  < 0 for  all  wavenumbers  a.  The  critical  Reynolds  number  Rc  is  that  for 
which  the  imaginary  part  of  5 first  vanishes  at  a critical  wavenumber  ac  as  R is  increased. 

The  use  of  the  Chebyshev  tau  method  to  discretize  the  equations  produces  two  spurious 
eigenvalues  with  Im(s)  > 0;  results  for  the  standard  test  case  a = 1 and  R = 10,000  [1]  are 
given  in  Table  II.  The  results  were  obtained  using  single  precision  on  the  CDC  Cyber  205  at 
NIST  with  the  IMSL  routine  EIGZC  [12]. 

The  spurious  eigenvalues  may  be  removed  by  employing  the  vorticity-streamfunction 
formulation  [2,7]  or  by  employing  a Galerkin  [2]  or  modified  Galerkin  approach  [6].  As  in 
the  model  problem,  Zebib’s  approaches  are  expected  to  lead  to  the  same  eigenvalues  as  the 
Chebyshev  tau  [4]  or  Chebyshev  Galerkin  [6]  methods  (in  the  absence  of  rounding  errors). 

As  another  alternative,  a direct  modification  of  the  usual  Chebyshev  tau  formulation  is 
suggested  by  the  results  of  section  2.2.  If  the  equations  are  written  in  the  form 

[ipzzzz  - 2 a2V>*2  + a4xp]/(—iaR)  + U(z)(^zz  - a2?/')  - Uzzip  = s(^zz  - a2il>), 
then  the  Chebyshev  tau  formulation  produces  N — 3 linear  equations 

[of  -2a2a(V  + a\n\l(-iaR)  + ...  = s(a[2)  - a2an)  n = 0, . . . , N - 4,  (19) 
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together  with  the  four  boundary  conditions  (17b),  where  for  simplicity  we  have  omitted  in 
(19)  the  convolution  sums  involving  U(z)  whose  form  is  unimportant  for  this  discussion.  To 
obtain  the  suggested  modification,  this  expression  is  approximated  instead  by  the  formula 

[a{4)  - 2 a2a^]  + a4an]/{-iaR ) + . . . = s{a(2)  - a2an)  n = 0, . . . , N - 4,, 

where  the  left  hand  side  is  unchanged  and  a ^ is  given  by  (8);  that  is,  the  last  two  columns 
in  the  matrix  that  represents  the  second  derivative  are  set  to  zero.  As  shown  in  Table  III, 
this  modification  also  serves  to  eliminate  the  two  spurious  eigenvalues  in  this  case  as  well, 
with  essentially  no  loss  of  accuracy  in  the  other  eigenvalues.  It  is  not  surprising  that  the 
modification  successfully  eliminates  the  spurious  modes,  since  the  Orr-Sommerfeld  equation 
differs  from  the  model  problem  only  in  the  lower-order  derivative  terms  on  either  side  of  the 
equation. 

Due  to  differing  treatment  of  these  lower-order  terms,  this  modification  of  the  usual 
tau  formulation  of  the  Orr-Sommerfeld  equation  is  not  equivalent  to  either  the  vorticity- 
streamfunction  formulation  or  the  Galerkin  formulation;  the  corresponding  eigenvalues  that 
are  computed  for  a given  N are  similar  but  not  identical.  Our  modified  tau  method  is  more 
efficient  than  the  vorticity-streamfunction  formulation  [2],  which  doubles  the  size  of  the  sys- 
tem, or  a Galerkin  formulation  [1,6],  which  is  more  awkward  to  implement.  Eliminating  half 
of  the  unknowns  to  reduce  the  size  of  the  system  in  the  vorticity-streamfunction  formulation 
[7]  also  produces  a more  complicated  system  of  equations. 

5.  Acknowledgements 

The  authors  are  grateful  for  helpful  discussions  with  S.  Coriell,  J.  Gary,  A.  Pearlstein,  and  A. 
Zebib.  The  first  author  (GBM)  acknowledges  partial  support  by  the  Microgravity  Sciences 
and  Applications  Program,  NASA.  The  second  author  (BTM)  was  supported  by  an  National 
Research  Council  Postdoctoral  Research  Fellowship. 


-12- 


References 

[1]  S.  A.  Orszag,  J.  Fluid  Mech.  50,  689  (1971). 

[2]  D.  Gottlieb  and  S.  A.  Orszag,  Numerical  A nalysis  of  Spectral  Methods:  Theory  and  Ap- 
plications (CBMS-NSF  Regional  Conference  Series  in  Applied  Mathematics  26,  SIAM, 
Philadelphia,  1977). 

[3]  J.  Gary  and  R.  Helgason,  J.  Comput.  Phys.  5,  169  (1970). 

[4]  A.  Zebib,  J.  Comput.  Phys.  53,  443  (1984). 

[5]  B.  Brenier,  B.  Roux,  and  P.  Bontoux,  J.  Theor.  Appl.  Mech.  5,  95  (1986). 

[6]  A.  Zebib,  J.  Comput.  Phys.  70,  521  (1987). 

[7]  D.  R.  Gardner,  S.  A.  Trogdon  and  R.  W.  Douglass,  J.  Comput.  Phys.  80,  137  (1989). 

[8]  B.  S.  Garbow,  J.  M.  Boyle,  J.  J.  Dongarra,  and  C.  B.  Moler,  Matrix  Eigensy.stem  Routine 
- EISPACK  Guide  Extension , (Springer-Verlag,  New  York,  1977). 

[9]  C.  B.  Moler  and  G.  W.  Stewart,  SIAM  J.  Numer.  Anal.  10,  241  (1973). 

[10]  G.  W.  Stewart,  in  Recent  Advances  in  Numerical  Analysis , 1978,  edited  by  C.  de  Boor 
and  G.  H.  Golub,  eds.,  (Academic  Press,  New  York, 1978),  p.  193. 

[11]  L.  M.  Delves  and  J.  L.  Mohamed,  Computational  Methods  for  Integral  Equations , (Cam- 
bridge University  Press,  New  York,  1985),  p.  212. 

[12]  T.  J.  Aird,  “ThelMSL  Library,”  in  Sources  and  Development  of  Mathematical  Software , 
edited  by  E.  Cowell,  (Prentice  Hall,  New  Jersey,  1984),  p.  264. 


-13- 


Table  I.  Eigenvalues  for  the  model  problem 


Chebyshev  tau 

[2] 

Present 

results 

N 

<7i 

<75 

<7  max 

<7l 

<75 

10 

-9.8696598 

-189.63800 

4,272. 

-9.8695970 

-97.95740 

15 

-9.8696044 

-89.54550 

29,439. 

-9.8696044 

-88.84327 

20 

-9.8696044 

-88.82644 

111,226. 

-9.8696044 

-88.82644 

25 

-9.8696044 

-88.82644 

294,697. 

-9.8696044 

-88.82644 

30 

-9.8696044 

-88.82644 

652,722. 

-9.8696044 

-88.82644 

35 

-9.8696044 

-88.82644 

1,255,298. 

-9.8696044 

-88.82644 

Exact 

-9.8696044 

-88.82644 

-9.8696044 

-88.82644 

Table  II.  Occurrence  of  spurious  eigenvalues  in  the  Chebyshev  tau  method.0 


No. 

Equation 

Boundary  Conditions 

Spurious 

1 

z/d2) 

= CTIp 

ip{±  1)  - 0 

None 

2 

l/d2) 

--  crip 

ipW{±l)  = 0 

None 

3 

- crip 

Ip{±  1)  = ^(1)(±1)  = 0 

None 

4 

1p^ 

--  tjipW 

ip{±  1)  - 0(1)(±1)  = 0 

Two 

5 

0(4) 

— (Tip d 

ip(dzl)  = ip^(±l)  = 0 

None6 

6 

d’(6) 

= crip 

V>(±1)  — V^(±l)  = 0(2)(d:l)  = 0 

None 

7 

0(6) 

— crip O 

lp(±l)  - ^^(±1)  = lp^(dc.l)  = 0 

Two 

8 

0(6) 

— crip d 

t/>(±l)  — ip^(±l)  = ip^2\±  1)  = 0 

Fourc 

a.  Here  we  denote  derivatives  of  ^ by  a superscript  in  parentheses. 

b.  May  be  reduced  to  a Dirichlet  problem  for  ip^2\ 

c.  Two  complex  conjugate  pairs  with  Re(cr)  < 0,  but  Im(a)  ^ 0. 
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Table  III.  First  four  eigenvalues  for  the  Orr-Sommerfeld  equation. 


a = 1 and  R — 

104. 

N 

Chebvshev  tau 

Present  results 

26 

0.08079254  + 35.0190231  i 
0.08896405  + 29.3934954  i 
0.23713751  + 0.00563644  i 
0.76774613  - 0.00334424  i 

0.23627968  + 0.00445813  i 
0.89244759  - 0.02021587  i 
0.96105877  - 0.02035856  i 
0.90085138  - 0.02192943  i 

38 

0.05429957  + 178.047201  i 
0.05782029  + 158.755478  i 
0.23752676  + 0.00373427  i 
0.96383565  - 0.03503110  i 

0.23752985  + 0.00373031  i 
0.96336687  - 0.03521228  i 
0.96270591  - 0.03619040  i 
0.90798372  - 0.04519238  i 

50 

0.04092607  + 563.551154  i 
0.04288092  + 517.476254  i 
0.23752648  + 0.00373967  i 
0.96462865  - 0.03516827  i 

0.23752648  + 0.00373967  i 
0.96462731  - 0.03516958  i 
0.96464022  - 0.03518657  i 
0.27720546  - 0.05089517  i 

From  [l] 

0.23752649  + 0.00373967  i 
0.96463092  - 0.03516728  i 
0.96464251  - 0.03518658  i 
0.27720434  - 0.05089873  i 

0.23752649  + 0.00373967  i 
0.96463092  - 0.03516728  i 
0.96464251  - 0.03518658  i 
0.27720434-  0.05089873  i 
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